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Abstract 

A combination  of  cartesian-grid  methods  and  volume-of- 
fluid  methods  is  used  to  simulate  breaking  waves  around 
ships  and  the  resulting  hydrodynamic  forces.  A surface 
panelization  of  a ship  hull  is  used  as  input  to  automat- 
ically generate  an  immersed-boundary  representation  of 
the  geometry  on  a cartesian  grid.  No  additional  gridding 
beyond  what  is  already  used  in  potential-flow  methods 
and  hydrostatics  calculations  is  required.  The  volume- 
of-fluid  portion  of  the  numerical  algorithm  is  used  to  cap- 
ture the  free-surface  interface,  including  the  breaking  of 
waves,  the  formation  of  spray,  and  the  entrainment  of  air. 
The  numerical  scheme  is  implemented  on  a parallel  com- 
puter. The  numerical  simulations  are  compared  to  analyt- 
ical solutions  and  experimental  measurements.  Together, 
the  ease  of  input  and  usage,  the  ability  to  model  and  re- 
solve complex  free-surface  phenomena,  and  the  speed  of 
the  numerical  algorithm  provide  a robust  capability  for 
simulating  the  free-surface  disturbances  near  a ship. 


Introduction 

The  Numerical  Flow  Analysis  (NFA)  code  provides 
turnkey  capabilities  to  model  breaking  waves  around 
a ship,  including  both  plunging  and  spilling  breaking 
waves,  the  formation  of  spray,  and  the  entrainment  of  air. 
The  basic  details  of  the  numerical  algorithm  are  provided 
in  Sussman  & Dommermuth  (2001),  and  Dommermuth 
et  al.  (2004),  and  Dommermuth  et  al.  (2006).  An  eval- 
uation of  NFAs  capabilities  is  provided  in  Wilson  et  al. 
(2006).  Here,  the  numerical  formulation  of  NFA  will  be 
discussed  in  detail,  including  the  treatment  of  the  hull 
boundary  condition  and  the  update  of  the  free-surface  el- 
evation. The  calculation  of  hydrodynamic  forces  and  the 
modeling  of  nonlinear  incident  waves  will  be  used  to  il- 
lustrate the  accuracy  and  stability  of  the  numerical  for- 


mulation. 

NFA  uses  a cartesian-grid  formulation  with 
immersed-body  and  volume-of-fluid  (VOF)  methods. 
The  governing  equations  are  formulated  on  a cartesian 
grid  thereby  eliminating  complications  associated  with 
body-fitted  grids.  The  sole  geometric  input  into  NFA 
is  a surface  panelization  of  the  ship  hull.  No  additional 
gridding  beyond  what  is  already  used  in  potential-flow 
methods  and  hydrostatics  calculations  is  required.  The 
ease  of  input  in  combination  with  a flow  solver  that  is 
implemented  using  parallel-computing  methods  permit 
the  rapid  turn  around  of  numerical  simulations  of 
complex  interactions  between  free  surfaces  and  ships. 

Based  on  Colella  et  al.  (1999)  and  Sussman  & Dom- 
mermuth (2001),  free-slip  boundary  conditions  are  im- 
posed on  the  surface  of  the  ship  hull.  A surface  repre- 
sentation of  the  ship  hull  is  used  as  input  to  construct 
a three-dimensional  signed-distance  representation.  The 
signed-distance  function  is  used  to  calculate  how  the  ship 
hull  cuts  through  the  cartesian  grid.  The  resulting  frac- 
tional areas  and  volumes  are  then  used  in  a finite-volume 
method  to  project  the  flow  velocity  onto  a solenoidal 
field  and  to  impose  the  hull  boundary  conditions.  The 
fractional  areas  and  volumes  of  very  small  cells  are 
merged  to  improve  the  conditioning  of  the  Poisson  solver 
(Mampaey  & Xu  1995).  Free-slip  boundary  conditions 
are  also  imposed  along  the  sides  and  at  the  top  and  bot- 
tom of  the  computational  domain.  At  the  exit,  Orlanski- 
like  boundary  conditions  are  imposed  (Orlanski  1976). 
At  the  entrance,  a wavemaker  is  used  to  generate  ambi- 
ent waves.  The  ambient  waves  are  accurate  to  second  or- 
der in  wave  steepness.  The  water-particle  velocity  in  the 
water  and  air,  and  the  free-surface  elevation  is  assigned 
at  the  entrance  to  the  computational  domain. 

The  grid  is  stretched  along  the  cartesian  axes  using 
one-dimensional  elliptic  equations  to  improve  resolution 
near  the  ship  hull  and  the  free-surface  interface.  Away 


from  the  ship  and  the  free-surface  interface,  where  the 
flow  is  less  complicated,  the  mesh  is  coarser.  Details  of 
the  grid  stretching  algorithm,  which  uses  weight  func- 
tions that  are  specified  in  physical  space,  are  provided  in 
Knupp  & Steinerg  (1993). 

The  VOF  portion  of  the  numerical  algorithm  is  used 
to  track  the  free-surface  interface,  including  the  large- 
scale  effects  of  breaking  waves,  spray  formation  and 
air  entrainment.  The  interface  tracking  of  the  free  sur- 
face is  second-order  accurate  in  space  and  time.  At 
each  time  step,  the  position  of  the  free  surface  is  re- 
constructed using  piece-wise  planar  surfaces  (Rider  et  al. 
1994,  Gueyffier  et  al.  1999).  The  advection  portion  of  the 
VOF  algorithm  uses  an  operator-split  method  (Puckett 
et  al.  1997).  The  advection  algorithm  implements  a cor- 
rection to  improve  mass  conservation  when  the  flow  is 
not  solenoidal  due  to  numerical  errors. 

The  convective  terms  in  the  momentum  equations 
are  treated  using  a slope-limited,  third-order  QUICK 
scheme  as  discussed  in  Leonard  (1997).  A Smagorin- 
sky  turbulence  model  is  also  implemented.  There  are  no 
special  treatments  required  to  model  either  the  flow  sepa- 
ration at  the  transom  or  the  wave  overturning  at  the  bow. 
A second-order,  variable-coefficient  Poisson  equation  is 
used  to  project  the  velocity  onto  a solenoidal  field.  A pre- 
conditioned conjugate-gradient  method  is  used  to  solve 
the  Poisson  equation. 

NFA  is  written  in  Fortran  90.  The  governing  equa- 
tions are  solved  using  a domain-decomposition  method. 
The  domains  are  distributed  over  the  nodes  of  a paral- 
lel computer.  Communication  between  processors  on  the 
Cray  XT3  is  performed  using  either  Crays  shared  mem- 
ory access  library  (SHMEM).  NFA  also  runs  on  clusters 
using  MPI.  The  CPU  requirements  are  linearly  propor- 
tional to  the  number  of  grid  points  and  inversely  propor- 
tional to  the  number  of  processors.  Together,  the  ease  of 
input  and  usage,  the  ability  to  model  and  resolve  complex 
free-surface  phenomena,  and  the  speed  of  the  numerical 
algorithm  provide  a robust  capability  for  simulating  the 
free-surface  disturbances  near  a ship. 


Formulation 


Mj  and  Xi  are  normalized  by  Ug  and  Lg,  which  denote 
the  free-stream  velocity  and  the  length  of  the  body,  re- 
spectively. The  frame  of  reference  is  fixed  with  respect 
to  the  initial  position  of  the  ship’s  fore  perpendicular.  A 
right-handed  coordinate  system  is  used  with  x positive 
forward  and  z positive  upward. 

Following  a procedure  that  is  similar  to  Rider  et  al. 
(1994),  we  let  cf)  denote  the  fraction  of  fluid  that  is  inside 
a cell.  By  definition,  ((>  = 0 for  a cell  that  is  totally  filled 
with  air,  and  <()  = 1 for  a cell  that  is  totally  filled  with 
water. 

The  advection  of  (/>  is  expressed  as  follows; 
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<5  is  a sub-grid-scale  flux  which  can  model  the  entrain- 
ment of  gas  into  the  liquid.  Dommermuth  et  al.  (1998) 
provide  details  of  a sub-grid  model  that  is  appropriate  for 
interface  capturing  methods  that  allow  mixing  of  air  and 
water.  Q = 0 for  the  present  formulation. 

Let  pe  and  pe  respectively  denote  the  density  and 
dynamic  viscosity  of  water.  Similarly,  pg  and  pg  are  the 
corresponding  properties  of  air.  The  flows  in  the  water 
and  the  air  are  governed  by  the  Navier-Stokes  equations: 
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where  Re  = peUgLo/pe  is  the  Reynolds  number  and 
= Ug/{gLo)  is  the  Froude  number,  g is  the  accel- 
eration of  gravity.  5,3  is  the  Kronecker  delta  function. 
P is  the  pressure.  As  described  in  Dommermuth  et  al. 
(1998),  Tij  is  the  subgrid-scale  stress  tensor.  S'y  is  the 
deformation  tensor: 
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p and  p are  respectively  the  dimensionless  variable  den- 
sities and  viscosities: 


Consider  turbulent  flow  at  the  interface  between  air  and 
water.  A two-phase  formulation  of  the  Navier-Stokes 
equations  is  used  to  model  a ship  moving  with  constant 
forward  speed  with  incident  waves.  Let  m,  denote  the 
three-dimensional  velocity  field  as  a function  of  space 
(xi  = (x,y,z))  and  time  (t).  Vi  is  the  velocity  of  the 
ship.  Vi  includes  the  effects  of  rigid-body  translation  and 
rotation.  Let  U,  denote  the  free-stream  current.  For  an 
incompressible  flow,  the  conservation  of  mass  gives 


p{4>)  = A + (l-A)H(^) 

= ??  + (1  - , (5) 

where  A = pg/ p(  and  77  = Pg/pi  are  the  density  and  vis- 
cosity ratios  between  air  and  water.  For  a sharp  interface, 
with  no  mixing  of  air  and  water,  H is  a step  function.  In 
practice,  a mollified  step  function  is  used  to  provide  a 
smooth  transition  between  air  and  water. 

A no-flux  condition  is  imposed  on  the  surface  of  the 
ship  hull: 


(1) 


UiTli  + UiTli  = ViTli 


(6) 


where  n*  denotes  the  normal  to  the  ship  hull  that  points 
into  the  fluid. 

As  discussed  in  Dommermuth  et  al.  (1998),  the  di- 
vergence of  the  momentum  equations  (3)  in  combina- 
tion with  the  conservation  of  mass  ( 1 ) provides  a Poisson 
equation  for  the  dynamic  pressure: 

axi  p axi 

where  S is  a source  term.  As  shown  in  the  next  sec- 
tion, the  pressure  is  used  to  project  the  velocity  onto  a 
solenoidal  field. 

Numerical  Time  Integration 


and  the  volume  fraction  is  advanced  to  complete  the  al- 
gorithm; 

(13) 

Gridding 

Along  the  cartesian  axes,  one-dimensional  stretch- 
ing is  performed  using  a differential  equation.  Let  x de- 
note the  position  of  the  grid  points  in  physical  space,  and 
let  ^ denote  the  position  of  the  grid  points  in  a mapped 
space.  As  shown  by  Knupp  & Steinerg  (1993),  the  dif- 
ferential equation  that  describes  grid  stretching  in  one  di- 
mension is  as  follows; 


Based  on  Sussman  (2003a),  a second-order  Runge- 
Kutta  scheme  is  used  to  integrate  with  respect  to  time  the 
field  equations  for  the  velocity  field.  Here,  we  illustrate 
how  a volume  of  fluid  formulation  is  used  to  advance 
the  volume-fraction  function.  Similar  examples  are  pro- 
vided by  Rider  et  al.  (1994).  During  the  first  stage  of  the 
Runge-Kutta  algorithm,  a Poisson  equation  for  the  pres- 
sure is  solved: 
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where  Ri  denotes  the  nonlinear  convective,  hydrostatie, 
viscous,  and  sub-grid-scale  terms  in  the  momentum 
equations,  uf,  and  P*  are  respectively  the  veloe- 
ity  components,  density,  and  pressure  at  time  step  k.  At 
is  the  time  step. 

For  the  next  step,  this  pressure  is  used  to  projeet  the 
velocity  onto  a solenoidal  field.  The  first  prediction  for 
the  velocity  field  (u*)  is 

The  volume  fraction  is  advanced  using  a volume  of  fluid 
operator  (VOF): 
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where  w{x)  is  a weight  function  that  is  specified  in  phys- 
ical space.  For  example,  suppose  the  grid  spacing  is  con- 
stant but  different  for  x < Xo  and  x > xi.  Between 
To  < X < xi,  there  is  a transition  zone  from  one  grid 
spacing  to  the  next.  Then  the  following  weight  function 
may  be  used  to  describe  this  distribution  of  grid  points: 
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Using  this  approach,  multiple  zones  of  grid  clustering 
may  be  specified.  For  example,  along  the  x-axis  (xi  in 
indical  notation),  grid  points  may  be  clustered  near  the 
bow  and  stern.  For  the  y-axis  (x2  in  indical  notation), 
grid  points  are  clustered  near  the  centerline  out  beyond 
the  half  beam.  Finally,  for  the  z-axis  (X3  in  indical  no- 
tation), grid  points  are  clustered  near  the  mean  waterline 
in  a region  that  is  between  the  top  and  bottom  of  the  ship 
hull.  Note  that  equation  14,  is  a nonlinear  equation  that 
is  solved  iteratively. 


,/,*=/-  VOF  At)  (10) 


Details  of  the  VOF  operator  are  provided  later.  A Pois- 
son equation  for  the  pressure  is  solved  again  during  the 
second  stage  of  the  Runge-Kutta  algorithm: 
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Ui  is  advanced  to  the  next  step  to  complete  one  cycle  of 
the  Runge-Kutta  algorithm: 
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Enforcement  of  Body  Boundary  Conditions 

A no-flux  boundary  condition  is  imposed  on  the  sur- 
face of  the  body  using  a finite-volume  technique.  A 
signed  distance  function  ^ is  used  to  represent  the  body. 
V'  is  positive  outside  the  body  and  negative  inside  the 
body.  The  magnitude  of  tp  is  the  minimal  distance  be- 
tween the  position  of  ip  and  the  surface  of  the  body,  ip  is 
zero  on  the  surface  of  the  body,  ip  is  calculated  using  a 
surface  panelization  of  the  hull  form.  Green’s  theorem  is 
used  to  indicate  whether  a point  is  inside  or  outside  the 
body,  and  then  the  shortest  distance  from  the  point  to  the 
surface  of  the  body  is  calculated.  Triangular  panels  are 


used  to  discretize  the  surface  of  the  body.  The  shortest 
distance  to  the  surface  of  the  body  can  occur  on  either 
a surface,  edge,  or  vertice  of  a triangular  panel.  Details 
associated  with  the  calculation  of  i/'  are  provided  in  Suss- 
man  & Dommermuth  (2001). 


Cells  near  the  ship  hull  may  have  an  irregular  shape, 
depending  on  how  the  surface  of  the  ship  hull  cuts  the 
cell.  On  these  irregular  boundaries,  the  finite-volume  ap- 
proach is  used  to  impose  free-slip  boundary  conditions. 
Let  Sh  denote  the  portion  of  the  cell  whose  surface  is  on 
the  body,  and  let  So  denote  the  other  bounding  surfaces 
of  the  cell  that  are  not  on  the  body.  Gauss’s  theorem  is 
applied  to  the  volume  integral  of  equation  8: 
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Here,  rii  denotes  the  components  of  the  unit  normal  on 
the  surfaces  that  bound  the  cell.  Based  on  equation  9,  a 
Neumann  condition  is  derived  for  the  pressure  on  Sb  as 
follows: 
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The  Neumann  condition  for  the  velocity  (6)  is  substituted 
into  the  preceding  equation  to  complete  the  Neumann 
condition  for  the  pressure  on  Sb'. 


dP^-  U*rii  v*rii 


p{4>’^')  dxi 


At 


At 


+ 


At 


+ RtUi  . (18) 


This  Neumann  condition  for  the  pressure  is  substituted 
into  the  integral  formulation  in  equation  16: 


ds 
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dxi 


This  equation  is  solved  using  the  method  of  fractional 
areas.  Details  associated  with  the  calculation  of  the 
area  fractions  are  provided  in  Sussman  & Dommermuth 
(2001)  along  with  additional  references.  Cells  whose  cut 
volume  is  less  than  25%  of  the  full  volume  of  the  cell  are 
merged  with  neighbors.  The  merging  occurs  along  the 
direction  of  the  steepest  gradient  of  the  signed-distance 
function  'ip.  This  improves  the  conditioning  of  the  Pois- 
son equation  for  the  pressure.  As  a result,  the  stability  of 
the  projection  operator  for  the  velocity  is  also  improved 
(see  equations  9 and  12). 


Interface  reconstruction  and  advection 


In  our  VOF  formulation,  the  free  surface  is  recon- 
structed from  the  volume  fractions  using  piece-wise  lin- 
ear polynomials.  The  reconstruction  is  based  on  algo- 
rithms that  are  described  by  Gueyffier  et  al.  (1999).  The 
surface  normals  are  estimated  using  weighted  central  dif- 
ferencing of  the  volume  fractions.  A similar  algorithm  is 
described  by  Pilliod  & Puckett  (1997).  Near  the  body, 
care  must  be  taken  to  use  cells  whose  volume  fraction  is 
exterior  to  the  body  in  the  calculation  of  the  normal  to 
the  free-surface  interface.  The  advection  portion  of  the 
algorithm  is  operator  split,  and  it  is  based  on  similar  al- 
gorithms reported  in  Puckett  et  al.  (1997).  Major  differ- 
ences between  the  present  algorithm  and  earlier  methods 
include  special  treatments  to  account  for  the  body  and  to 
alleviate  mass-conservation  errors  due  to  the  presence  of 
non-solenoidal  velocity  fields. 

Let  Fi  denote  flux  through  the  faces  of  a cell.  Fi  is 
expressed  in  terms  of  the  relative  velocity  (uj  -\-Ui-Vi) 
and  the  areas  of  the  faces  of  the  cell  (Aj)  that  are  cut  by 
the  ship  hull: 

Fi  = Ai  (ui  + Ui-  Vi)  . (20) 


If  the  ship  hull  does  not  cut  the  cell,  then  A*  correspond 
to  the  surface  areas  that  bound  the  cell.  Near  the  ship 
hull,  Ai  is  some  fraction  of  the  surface  areas  that  bound 
the  cell.  Note  that  A*  = 0 inside  the  ship  hull.  Based  on 
an  application  of  Gauss’s  theorem  to  the  volume  integral 
of  Equation  1 and  making  use  of  Equation  6: 

F+-F-=0,  (21) 

where  F^  is  the  flux  on  the  positive  i-th  face  of  the  cell 
and  F~  is  the  flux  on  the  negative  i-th  face  of  the  cell. 
Due  to  numerical  errors,  equation  21  is  not  necessarily 
satisfied.  Let  € denote  the  resulting  numerical  error  for 
any  given  cell.  For  each  cell  whose  flux  is  not  conserved, 
a correction  is  applied  prior  to  performing  the  VOF  ad- 
vection. For  example,  the  following  reassignment  of  the 
flux  along  the  vertical  direction  ensures  that  the  redefined 
flux  is  conserved: 
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As  the  ship  moves  due  to  force  motions  or  in  re- 
sponse to  incident  waves,  the  signed-distance  function 
V'  and  the  fractional  areas  and  volumes  are  recalculated 
each  time  step.  In  the  case  of  forced  motions,  the  signed- 
distance  function  can  be  precalculated  over  the  range  of 
motion,  then  interpolation  can  be  used  during  the  simula- 
tion to  calculate  the  signed-distance  function  to  improve 
efficiency. 


Based  on  this  new  flux,  new  relative  velocities  are  de- 
fined on  the  faces  of  the  cell: 


AxjPj 
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(23) 


where  Axi  is  the  grid  spacing  and  V = AxiAx2Axs 
is  the  volume  of  the  cell.  Away  from  the  ship  hull,  rii 


is  the  relative  velocity  plus  a corrective  term  to  conserve 
mass.  Inside  the  ship  hull,  Ui  = 0 because  A*  = 0.  Near 
the  ship  hull,  Ui  is  scaled  by  the  fraction  of  area  that  is 
cut  by  the  presence  of  the  ship  hull,  u-i  is  continuous 
across  the  faces  of  the  cells  along  x-  and  y-axes,  but 
discontinuous  across  the  faces  along  the  z-axis  because 
in  this  particular  example  that  is  the  axis  where  the  flux 
has  been  corrected. 

Equation  2 is  operator  split.  A dilation  term  is  added 
to  ensure  that  the  volume  fraction  remains  between  0 < 
(f)  < 1 during  each  stage  of  the  splitting  (Puckett  et  al. 
1997).  The  resulting  discrete  set  of  equations  for  the  first 
stage  of  the  time-stepping  procedure  is  provided  below: 
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Ti  denotes  VOF  advection  based  on  the  uncut  areas  of 
the  faces  of  the  cell.  As  an  example,  for  a cell  that  is  full 
of  water,  ,Fi(ui,  (/),  Af)  = <fmiAtAx2Ax3.  The  dila- 
tion term  is  treated  explicitly  in  the  first  two  parts  of  the 
operator-split  algorithm  and  implicitly  in  the  last  part  of 
the  preceding  equation.  Note  that  the  order  of  the  split- 
ting is  alternated  from  time  step  to  time  step  to  preserve 
second-order  accuracy. 


Interface  Visualization 

The  free-surface  interface  that  is  reconstructed  from 
the  volume  fractions  is  most  often  calculated  and  visual- 
ized using  commercial  codes.  Specifically,  commercial 
codes  calculate  the  0.5  isosurface  of  the  volume -fraction 
function  </>.  The  free-surface  interface  that  is  calculated 
from  the  0.5  isosurface  is  different  from  the  free-surface 
interface  that  is  reconstructed  from  the  volume  fractions. 
To  illustrate  this  point,  consider  a cell  whose  volume 
fraction  (j)o  is  between  half  full  and  full,  0.5  < </>o  < 1. 
Let  Az  denote  the  height  of  the  cell.  Assume  that  the 
free-surface  interface  is  horizontal  and  that  all  the  fluid 
is  sitting  in  the  bottom  of  the  cell.  Then  the  height  of  the 
free-surface  interface  above  the  bottom  of  the  cell  based 


on  VOF  reconstruction  is  as  follows: 


Tj  = (t>oAz  . (25) 

In  contrast,  if  the  cell  above  is  filled  with  air.  then  based 
on  the  0.5  isosurface,  the  height  of  the  free-surface  inter- 
face is 


The  maximum  difference  between  equations  25  and  26 
occurs  when  = 3/4.  The  error  at  this  point  is  about 
1 1 % higher  for  0.5  isosurface  relative  to  VOF  reconstruc- 
tion. If  the  volume  fraction  is  less  than  < 0.5  within 
a cell  and  its  neighbors,  then  the  0.5  isosurface  does  not 
even  exist.  This  is  problematic  in  visualizations  of  turbu- 
lent flows  with  lots  of  spray  because  droplets  and  sheets 
of  water  can  suddenly  appear  and  disappear. 


Radiation  Conditions 


Exit  boundary  conditions  are  required  in  order  to 
conserve  mass  and  flux.  For  ships  with  forward  speed,  an 
Orlanski-like  formulation  (Orlanski  1976)  provides  the 
necessary  radiation  condition. 


dui  ^ 

dt  ° dx 


= 0 


(27) 


Uo  is  the  free-stream  current  along  the  x-axis,  and 
is  the  water-particle  velocity  along  the  x-axis.  For  the 
other  components  of  velocity  and  the  volume  fraction, 
zero  gradients  are  imposed  at  the  exit  of  the  computa- 
tional domain: 


du2  du3  d(j) 

dx  dx  dx 


(28) 


Neumman  conditions  are  specified  for  the  pressure  in  a 
manner  that  is  very  similar  to  the  imposition  of  free-slip 
conditions  on  the  ship  hull  (see  equations  16  thru  19). 
Based  on  the  x-component  of  momentum. 
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Upon  substitution  of  equation  27  into  the  preceding 
equation,  the  following  Neumann  condition  is  derived  for 
the  pressure  at  the  exit  of  the  computational  domain: 
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This  equation  is  substituted  into  the  set  of  finite-volume 
equations  that  govern  the  pressure  (see  equation  19). 

Equations  27  thru  30  prevent  the  reflection  of  distur- 
bances back  into  the  interior  of  the  computational  do- 
main. However,  these  equations  do  not  guarantee  the 


conservation  of  mass.  In  order  to  conserve  mass,  a re- 
gridding  procedure  is  introduced  when  there  are  no  in- 
cident waves.  The  initial  volume  fraction  is  integrated 
for  the  grid  cells  that  are  on  the  leading  edge  of  the 
computational  domain.  This  integrated  quantity  is  used 
to  maintain  a constant  mean  water  level  at  the  entrance 
to  the  computational  domain.  At  the  end  of  each  time 
step,  changes  in  the  integrated  volume  fraction  are  cal- 
culated. Any  changes  in  the  integrated  volume  fraction 
are  eliminated  by  imposing  a vertical  velocity  that  brings 
the  mean  water  level  at  the  leading  edge  back  into  align- 
ment. The  velocity  correction  is  used  to  move  the  volume 
fractions  over  the  entire  computational  domain  either  up 
or  down,  depending  on  the  situation.  A VOF  method  is 
used  to  move  the  volume  fractions.  The  VOF  method  en- 
sures that  the  free-surface  interface  remains  sharp  during 
the  regridding  process. 

INITIAL  Transients 

Initial  transients  are  minimized  using  an  adjustment 
procedure.  An  analysis  of  adjustment  procedures  as  it 
applies  to  free-surface  problems  is  provided  in  Dommer- 
muth  (1994)  and  Dommermuth  (2000).  Let  f(t)  denote 
the  adjustment  factor  as  a function  of  time,  then  f{t)  and 
its  derivative  f'{t)  are  by  definition 

fit)  = 1 - exp(-(;^)2) 

fit)  = 2^exp(-(^)2)  , (31) 

where  To  is  the  adjustment  time. 

For  a ship  hull  that  is  oscillating  up  and  down,  the 
vertical  motion  {z)  and  vertical  velocity  {w)  of  the  ship 
hull  are 

z = Asiniujt)f{t) 

w = Auj  cos(ujt)  f (t)  + Asin(ujt)  f (t)  . (32) 

A is  the  amplitude  and  u>  is  the  frequency  of  the  ver- 
tical motion.  A similar  procedure  is  used  for  forced  roll 
and  pitch  motions.  The  free-stream  current  is  also  slowly 
ramped  up  to  speed. 

Enforcement  of  Courant  Conditions 

The  momentum  equations  are  integrated  in  time  us- 
ing an  explicit  Runge-Kutta  algorithm.  As  a result,  a 
Courant  condition  must  be  enforced  for  the  maximum 
total  velocity: 

/\  'j'  ■ 

\ui  -\-Ui\  < C (33) 

C is  a coefficient  that  ensures  that  the  Courant  condi- 
tion is  satisfied  for  both  the  momentum  equations  and 


the  VOF  advection.  Typically,  C = 0.45  in  the  numeri- 
cal results  that  are  presented  in  this  paper.  If  the  Courant 
condition  is  exceeded,  the  magnitude  of  the  velocity  is 
reduced  such  that  the  Courant  condition  is  satisfied.  This 
clipping  of  the  velocity  field  tends  to  occur  in  regions 
where  fine  spray  is  formed,  especially  in  the  rooster-tail 
region. 

Treatment  of  convective  terms 

The  convective  terms  in  the  momentum  equations 
(see  Equation  3)  are  calculated  using  a slope-limited, 
QUICK,  finite-difference  scheme  (Leonard  1997).  Spe- 
cial treatments  are  required  near  the  ship  hull.  One  pos- 
sibility is  to  use  one-sided  differencing.  However,  one- 
sided differencing  is  often  unstable.  Another  possibility 
is  to  extend  the  velocity  of  the  fluid  into  the  ship  hull. 
In  this  case,  setting  the  velocity  equal  to  zero  inside  the 
body  is  stable,  but  too  “sticky.”  Another  possibility  is  to 
extend  the  fluid  velocity  into  the  ship  hull  in  such  a man- 
ner that  the  no-flux  condition  is  met  right  at  the  ship  hull. 
The  interior  flow  that  meets  this  condition  is  as  follows: 

= [(^i  - tJj)nj]  Hi  , (34) 

where  recall  that  Vj  is  velocity  of  the  body,  t/j  is  the 
velocity  of  the  free-stream  current,  and  n,-  is  the  unit 
normal  that  points  along  gradient  of  the  signed-distance 
function  (i/>).  At  the  ship  hull,  UiUi  = vini  — UiTii  using 
this  formulation  of  the  interior  flow. 

Density  Smoothing 

The  density  as  a function  of  the  volume  fraction  is 
smoothed  using  a three-point  stencil  (1/4, 1/2, 1/4)  that  is 
applied  consecutively  along  each  of  the  cartesian  axes. 
This  improves  the  conditioning  of  the  Poisson  equation 
(Equations  8 & 1 1).  If  the  density  is  smoothed,  then  the 
same  smoothed  density  must  be  used  in  the  projection 
steps  (Equations  9 & 12). 

Calculation  of  forces 

The  forces  T)  acting  on  the  ship  hull  are  calculated 
in  two  parts  based  on  integration  of  the  normal  and  tan- 
gential stresses  over  the  surface  of  the  ship. 

Fi  = J dsPn-i  -f  J dsTi  , (35) 

where  the  normal  stress  is  expressed  in  terms  of  the  pres- 
sure P,  and  the  tangential  stress  n is  expressed  in  terms 
of  the  water-particle  velocity  just  outside  the  boundary 
layer  of  the  ship. 

Ti  = Cf)^p\U\(ui + Ui  - Vi)  , (36) 


where  cj  is  the  friction  coefficient  based  on  the  ITTC 
line. 

Cf  = 0.075(log(i?e)  - 2)-2  . (37) 

\U\  is  the  magnitude  of  the  water-particle  velocity. 

\U\  = \J {ui  + Ui-  Vi){ui  + Ui-  Vi)  . 

Note  that  {ui  + Ui  — Vi)rii  = 0.  Similarly,  the  moments 
Mi  are  expressed  as  follows: 

Mi  = j dsPeijkXjrik  -f  J dsCijtrjTk  , (38) 

where  eijk  is  the  Levi-Civita  function.  Vj  is  the  moment 
arm  and  eij^r^rik  is  the  indicial-notation  representaiton 
of  the  cross  product  of  the  moment  arm  with  the  unit 
normal. 

The  surface  integrals  for  the  forces  and  moments 
are  calculated  using  the  panelized  geometry.  The  den- 
sity, pressure,  and  water-particle  velocities  are  interpo- 
lated onto  the  panels.  The  integrands  are  evaluated  and 
summed  over  the  surface  of  the  body. 


Wavemaker 


Here,  we  highlight  the  formulation  of  a wavemaker. 
Let  r/(x,  t)  denote  the  free-surface  elevation  as  function 
of  position  x and  time  t,  then 

rj{x,t)  = af{t)cos{kx  + at) 

+ ^kaf{t)  cos(2fcx  + 2(jt)  , (39) 

where  a is  the  wave  amplitude,  k is  the  wavenumber,  and 
a is  the  encounter  frequency.  The  preceding  formula  is 
accurate  to  second  order  in  wave  steepness.  f(t)  is  an  ad- 
justment factor  that  slowly  ramps  up  the  wave  amplitude 
(see  Equation  3 1 ).  The  encounter  frequency  is  a function 
of  the  intrinsic  wave  frequency  w,  the  wavenumber,  and 
the  speed  of  the  free-stream  current  Uo  : 

a = uj  — kUo{t)  , (40) 

where  the  normalized  intrinsic  wave  frequency  is 

k 

ijp' = — irtsMaikd)  , (41) 

Ft 

where  d is  the  water  depth.  The  speed  of  the  ship  is 
slowly  ramped  up  from  rest  using  the  adjustment  factor. 

Uo{t)  = -fit)  ■ (42) 


For  z < rj , the  horizontal  and  vertical  components  of  the 
water-particle  velocity  are 

/ \ , cosh(A:(z -f  d))  „ , 

u[x,z,t)  = —aojj(t) ; — cos(kx  + at) 

smh(kd) 

. . , sinh(A:(2: -f  d))  . ,,  , 

w{x,z,t)  = —aojj{t) . sm[kx  + at)  . 
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(43) 


For  z > rj , the  horizontal  and  vertical  components  of  the 
air-particle  velocity  are 


u{x,  z,  t) 
w{x,  z,  t) 


aujf{t) 

aujf(t) 


cosh{k{z  — h)) 
smh(kh) 
sinh(A:(z  — h)) 
smh{kh) 


cos{kx  + at) 
sin(A:x  + at) 


(44) 


where  h is  the  height  of  the  air  above  the  free  surface. 
Using  this  formulation,  the  horizontal  water-particle  ve- 
locity is  discontinuous  across  the  air-water  interface,  and 
the  vertical  water-particle  velocity  is  continuous  across 
the  air-water  interface. 


Results 


5365  geometry  (Athena) 

A three-dimensional  numerical  simulation  that  uses 
850x192x128=  20,889,600  grid  points,  5x8x4=160  sub- 
domains,  and  1 60  nodes  has  been  performed  on  a Cray 
XT3.  The  length,  width,  depth,  and  height  of  the  com- 
putational domain  are  respectively  4.0,  1 .0,  1 .0,  0.5  ship 
lengths  (Lo).  Grid  stretching  is  employed  in  all  direc- 
tions. The  smallest  grid  spacing  is  0.0020  near  the  ship 
and  mean  waterline,  and  the  largest  grid  spacing  is  0.020 
in  the  far  field.  The  Froude  number  is  Fr  = 0.4316. 
Two  incident  wavelengths  are  considered:  A = 1/2  and 
A = 2.  In  both  cases,  the  wave  steepness  is  H/X  = 0.06, 
where  H = 2a  the  wave  height.  The  equations  for 
the  wavemaker  are  imposed  ahead  of  the  ship  over  the 
range  1 < x < 1.5.  Initial  transients  are  minimized  by 
slowly  ramping  up  the  free-stream  current  and  the  inci- 
dent wave  amplitude.  The  period  of  adjustment  is  To  = 
0.5.  For  this  simulation,  the  non-dimensional  time  step  is 
At=0.0005.  The  numerical  simulation  runs  10,000  time 
steps  corresponding  to  5 ship  lengths.  Each  simulation 
requires  about  80  hours  of  wall-clock  time.  The  numer- 
ical results  are  compared  to  experimental  measurements 
that  had  been  performed  at  the  Naval  Surface  Warfare 
Center  (Ratcliffe  & Dommermuth  2007). 

Figures  1 show  perspective  views  of  the  predicted 
free-surface  elevations  at  time  f = 5 for  the  two  cases. 
The  color  bar  denotes  free-surface  elevations  between  - 
0.02  to  0.02  for  the  short-wave  case  and  -0.04  to  0.04 
for  the  long-wave  case.  The  incident  wave  propagates 
from  the  inlet  to  the  exit  with  no  attenuation  in  amplitude. 
This  indicates  that  the  effects  of  numerical  dissipation 
are  minimal. 

The  predicted  and  measured  free-surface  elevations 
as  a function  of  time  at  are  shown  in  Figures  2.  The  pre- 
dicted free-surface  elevations  are  ramped  up  to  their  full 
height.  This  minimizes  transients  associated  with  start- 
ing up  the  numerical  wavemaker.  The  measured  free- 


surface  elevations  show  slight  irregularities  that  are  as- 
sociated with  limitations  with  the  wavemaker  at  DTMB. 

Figures  3 show  the  predicted  vertical  force  com- 
pared to  measurements.  The  displacement  has  been  sub- 
tracted out  from  the  results.  As  the  model  ramps  up  to 
full  speed,  a mean  suction  force  is  induced  on  the  model. 
The  oscillatory  portion  of  the  force  is  dominated  by  hy- 
drostatics. For  the  short-wave  case,  the  numerical  pre- 
dictions show  slightly  more  sinkage  than  had  been  mea- 
sured. The  oscillatory  portions  of  the  forces  agree  well. 

Figures  4 show  the  predicted  drag  force  compared  to 
measurements.  In  the  case  of  the  numerical  simulations, 
the  drag  is  initially  zero  because  the  model  is  ramped  up 
to  full  speed  from  zero  forward  speed.  The  primary  har- 
monics that  are  evident  in  the  plots  are  due  to  the  incident 
wave  forces.  We  speculate  that  the  higher  harmonics  that 
are  evident  in  the  laboratory  results,  which  are  especially 
evident  for  the  longer  wave,  are  due  to  vibrations  in  the 
structure  that  is  used  to  restrain  the  model.  Another  pos- 
sibility is  the  variation  in  the  incident  wave  amplitude  in 
the  experiments.  Structural  vibrations  could  be  assessed 
by  installing  accelerometers  on  the  model.  The  variation 
in  the  incident  wave  amplitude  could  be  accounted  for  in 
the  numerical  simulations  by  using  wave-probe  data  as 
input.  In  general,  numerical  predictions  and  laboratory 
measurements  agree  well. 

5514  and  5613  geometries 

Forced-motion  studies  of  models  5514  and  5613 
have  been  performed.  The  models  are  forced  to  move 
in  heave,  pitch,  and  roll  while  moving  with  constant 
forward  speed.  The  Froude  number  for  both  hull 
forms  is  Fr  = 0.3.  For  the  heave  and  pitch  stud- 
ies, a plane  of  symmetry  is  used  on  the  centerplane 
of  the  computational  domain.  The  heave  and  pitch 
simulations  use  680x192x128=16,711,680  grid  points 
and  4x8x4=128  sub-domains.  The  forced  roll  simu- 
lations use  680x384x128=33,423,360  grid  points  and 
4x8x4=128  sub-domains.  The  length,  width,  depth,  and 
height  of  the  computational  domain  are  respectively  3.0, 
1.0,  1.0,  0.5  ship  lengths  (Lo)  for  the  heave  and  pitch 
studies.  The  width  of  the  domain  is  doubled  for  the 
roll  simulations.  All  of  the  simulations  use  128  pro- 
cessors on  a Cray  XT3.  For  model  5514,  the  full-scale 
length  and  forward  speed  are  respectively  Lg  = 142.0m 
and  Uo  = 11.20m/s.  For  model  5613,  Lg  = 154.0m 
and  Ug  = 11.66m/s.  Table  1 summarizes  the  non- 
dimensional  amplitudes  and  frequencies  of  motion.  The 
heave  amplitudes  correspond  to  ranges  of  motion  that  are 
80%  of  the  mean  drafts.  Table  2 summarizes  the  time 
step  and  the  number  of  time  steps.  Smaller  time  steps 
are  used  for  the  roll  simulations  and  one  heave  case  due 
to  Courant-stability  considerations.  Several  animations 
of  these  simulations  have  been  prepared  at  the  flow  visu- 
alization center  at  ERDC.  The  animations  are  accessible 


1 5514  1 

1 5613  1 

amplitude 

frequency 

amplitude 

frequency 

Heave 

0.3667 

13.95 

0.02857 

14.53 

Pitch 

5° 

13.95 

5° 

14.53 

Roll 

65° 

8.525 

65° 

8.877 

Table  1 ; Forced-motion  parameters. 


1 5514  1 

1 5613  1 

At 

N 

At 

N 

Heave 

0.0005 

4,000 

0.00025 

10,000 

Pitch 

0.0005 

4,000 

0.0005 

4,000 

Roll 

0.00025 

8,000 

0.00025 

8,000 

Table  2:  Time  step  and  number  of  time  steps. 


at  http://www.saic.com/nfa/.  This  website  also  includes 
animations  of  the  Athena  moving  with  constant  forward 
speed  into  head  seas. 


Conclusions 

In  terms  of  progress,  it  is  interesting  to  consider  the  re- 
sults of  research  reported  in  earlier  ONR  symposiums. 
Dommermuth  et  al . ( 1 998)  study  the  flow  near  the  bow  of 
model  5415  using  a variable-density,  cartesian-grid  for- 
mulation. A body  force  is  used  by  Dommermuth  et  al. 
(1998)  to  impose  the  body  boundary  condition.  The  nu- 
merical results  of  Dommermuth  et  al.  (1998)  barely  cap- 
ture the  initial  onset  of  wave  overturning  near  the  bow. 
Sussman  & Dommermuth  (2001)  continue  to  develop 
interface  capturing  methods.  Once  again,  comparisons 
are  shown  to  the  bow  flow  of  model  5415.  The  results 
do  not  show  significant  improvement  over  their  earlier 
results.  However,  their  calculations  of  the  breakup  of 
a turbulent  spray  sheet  illustrate  a novel  application  of 
interface-capturing  methods.  Dommermuth  et  al.  (2004) 
use  two  methods  to  study  the  flow  around  model  5415,  a 
vertical  strut,  and  a bluff  wedge.  The  first  method  uses 
iree-slip  conditions  on  the  hull  in  combination  with  a hy- 
brid level-set  and  VOF  interface-capturing  method.  In 
addition,  adaptive  mesh  refinement  (AMR)  is  used  to  im- 
prove grid  resolution  near  the  hull  and  free-surface  inter- 
face. Their  preliminary  results  illustrate  the  efficiency 
of  AMR.  The  second  method  uses  body-force  and  VOF 
formulations  on  a cartesian  grid  with  no  grid  stretching. 
The  results  show  more  fine-scale  detail  than  the  earlier 
studies.  The  predicted  free-surface  elevations  compare 
well  with  experiments,  but  the  body-force  method  is  too 
“sticky”  because  too  much  fluid  is  dragged  with  the  ship 
hull.  Based  on  these  results,  Dommermuth  et  al.  (2006) 
use  free-slip  boundary  conditions  to  impose  the  body 
boundary  condition  to  reduce  stickiness.  The  VOF  algo- 
rithm is  generalized  to  include  free-slip  conditions  on  the 
ship  hull.  The  grid  is  stretched  along  the  cartesian  axes  to 


improve  grid  resolution.  Numerical  predictions  compare 
well  with  laboratory  measurements  of  ship  models  mov- 
ing with  constant  forward  speed.  The  present  research 
implements  forced  motions  and  ambient  waves.  As  the 
ship  moves,  the  portion  of  the  cartesian  grid  that  is  cut 
by  the  ship  hull  is  continuously  calculated.  The  inherent 
simplicity  of  the  cartesian-grid  formulation  is  retained 
even  though  the  ship  moves  through  the  grid.  Ambient 
waves  are  generated  by  forcing  the  flow  in  the  water  and 
the  air  based  on  analytic  solutions  for  two-phase  flows. 
The  procedure  is  illustrated  using  regular  waves,  but  it  is 
general  enough  to  model  a seaway. 
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Figure  3:  Model  5365  (Athena)  vertical  forces. 


Figure  4:  Model  5365  (Athena)  drag  forces. 


